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Abstract 

Immersed boundary methods for computing confined fluid and plasma 
flows in complex geometries are reviewed. The mathematical principle 
of the volume penalization technique is described and simple examples 
for imposing Dirichlet and Neumann boundary conditions in one dimen¬ 
sion are given. Applications for fluid and plasma turbulence in two and 
three space dimensions illustrate the applicability and the efficiency of the 
method in computing flows in complex geometries, for example in toroidal 
geometries with asymmetric poloidal cross-sections. 


1 Introduction 

Immersed boundary techniques, including penalization approaches, are nowa¬ 
days commonly employed to solve boundary or initial boundary value problems 
in complex geometries. They consist of embedding the original, possibly com¬ 
plex spatial domain inside a bigger domain having a simpler geometry, for exam¬ 
ple a Cartesian geometry, while keeping the boundary conditions approximately 
enforced thanks to new terms that are added to the equations. Historically, 
the penalization technique can be traced back to Courant :8] in the context of 
constrained optimization. Later, Saulyev |37] applied it in the context of fic¬ 
titious domain methods for immersed boundaries. In a classical paper 132] an 
immersed boundary method was proposed for computing the blood flow in a 
beating heart. Some history on immersed boundary techniques can be found in 
[3], and for detailed reviews we refer to the classical review papers |33] and (22j . 

In the current work we focus on one particular example, the volume penal¬ 
ization method [I] which, inspired by the physical intuition that a solid wall is 
similar to a vanishingly porous medium, uses the Brinkman-Darcy drag force 
as penalization term. The main advantage of such penalized equations is that 
they can be discretized independently of the geometry of the original problem, 
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since the latter has been encoded into the penalization terms. Such a simplifi¬ 
cation permits a massive reduction in solver development time, since it avoids 
the issues associated with the design and management of the grid, allowing for 
example the use of simple fast Fourier transform (FFT)-based spectral solvers in 
Cartesian geometries. Using spectral solvers no linear systems have to be solved, 
e.g., to impose the incompressibility of the magnetic and velocity fields. Fur¬ 
thermore the penalized equations are solved without introducing additional dis¬ 
cretization errors, and thus effects of numerical diffusion and dispersion present 
in low-order numerical schemes are avoided. Moreover, parallel FFT libraries 
(e.g., P3DFFT) are available on the current supercomputers. The advantage 
of the penalization method becomes even more substantial when the geometry 
is time-dependent, as in the case of moving obstacles, or when fluid-structure 
interaction is taken into account. 

In the following we give a non exhaustive overview of the development for 
different applications using volume penalization. Two-dimensional hydrody¬ 
namic flows in staggered or inline tube bundles have been computed in [58] and 
a dipole-wall collision benchmark using different numerical discretizations of 
the penalized Navier-Stokes equations has been obtained in M- Final states in 
two-dimensional decaying hydrodynamic turbulence in different geometries have 
been studied in [40] and flows past flat plates focusing on the influence of the 
plate’s geometry in [32] ■ Nguyen van yen et al. (2011) (22) computed dipole-wall 
collisions and derived the scaling of the energy dissipation in the large Reynolds 
number limit. They showed that it converges indeed to a finite value. The penal¬ 
ized wave equation has been analysed in m- Applications to moving obstacles 
and fluid-structure interaction problems can be found in [mnnim and [9]. In 
m a numerical study using the volume penalization method for impeller-driven 
von Karman flows has been performed. Different simulations of compressible 
flows in complex geometries using the volume penalization have been described 
in mm- An adaptive penalty method for incompressible Navier-Stokes has 
been proposed by Shirokoff and Nave [43] . 

The extension of penalization techniques to magnetohydrodynamic (MHD) 
flows is more recent. Simulations of the self-organization of confined plasma 
in two space dimensions studying the influence of the confinement geometry 
have been presented in HH3. The numerical method has been extended to 
three space dimensions and is described and benchmarked in detail in )25| . 
The spontaneous spin-up of plasma in toroidal geometries using the viscous- 
resistive MHD equations has been discovered in Morales et al. [21] ■ The effect 
of toroidicity in reversed field pinch devices is studied in [22]. Roberts et al. [35] 
investigated helically forced MHD flows in cylindrical geometries and showed 
the persistence of helical modes even when the dynamics becomes turbulent A 
penalty method for hyperbolic partial differential equations modeling the edge 
plasma transport in a tokamak has been proposed in [2j . 

The aim of this paper is to review the volume penalization method and to 
illustrate its potential for applications in fluid and plasma turbulence in complex 
domains. Typically Dirichlet boundary conditions are considered corresponding 
to imposing the value of the solution at the boundary. Extensions for deal¬ 
ing with Neumann boundary conditions, i.e., imposing values of the derivative 
of the solution at the boundary, have been proposed in mm based on pre¬ 
vious work in 33] . An alternative approach for imposing either Dirichlet or 
Neumann boundary conditions which is based on sharp interface methods has 
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been introduced in |25| . Considering simple examples in one space dimension 
allows an understanding and analysis of the convergence behavior of the pe¬ 
nalization techniques when the penalization parameter tends to zero. We show 
that, for a given numerical discretization of the penalized equations, there ex¬ 
ists a value of the penalization parameter, corresponding to a balance between 
penalization and discretization errors, below which no further gain in precision 
is achieved. These results shed light on the behavior of volume penalization 
schemes when solving the penalized equations, outline the limitations of the 
method, and give indications on how to choose the penalization parameter at 
least in simple test cases. Nevertheless for practical applications a series of com¬ 
putations may be necessary to determine the actual value of the parameters, as 
discussed, e.g., in HU. Finally, different illustrations will be given for hydro- 
and magnetohydrodynamic problems in the turbulent regime including a study 
of the spatio-temporal self-organization of visco-resistive magnetohydrodynam¬ 
ics in toroidal geometries with different poloidal cross-sections while imposing 
curl-free toroidal magnetic and electric fields. 

The outline of the paper is the following. First, we present a short primer 
on penalization. Then we describe in some detail the volume penalization to 
impose either Dirichlet or Neumann boundary conditions and we present some 
simple one-dimensional examples. Applications to fluid and plasma turbulence 
in two and three dimensions illustrate the properties of the method. Finally, we 
set out the conclusions together with some perspectives for future work. 


2 A short primer on penalization 


To illustrate the idea of volume penalization we consider a boundary (BVP) or 
initial boundary value problem (IBVP) for a partial differential equation in a 
domain flf C R d , written in the following abstract operator equation 

Lu = f for xGflf (1) 


completed with boundary conditions Bu = g at dflf and additional initial 
conditions in case of an IBVP. The differential operator L stands, for example, 
for the Laplace operator, L = —V 2 , or the Navier-Stokes or Maxwell operator. 
The boundary conditions can be of Dirichlet or Neumann type, e.g., u = g or 
du/dn = g , respectively, where n is the outer normal of the domain. Solving 
eq. Q numerically requires domain-fitted grids using, e.g., finite elements. The 
grid generation and the numerical solution of the resulting discretized problem 
can be demanding, see for example { 121 . 

An alternative are penalization methods which embed the problem posed 
in the complexly shaped domain f If into a larger simple domain, typically of 
rectangular shape f2, i.e., Clf C fl. The advantage is that fast solvers are 
available for such problems, which can furthermore be easily parallelized. 

In the Dirichlet case the penalized problem thus reads 


Lu v = f 


~x(Bu 

V 


g) for xGfl = D/Un s , 


( 2 ) 


where the boundary conditions have been included into the equation and an 
additional parameter, the small penalization parameter 77 > 0 has been intro¬ 
duced. All information about the geometry of Qf has been encoded into the 
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mask function \i which is defined as 


for x £ flf 
elsewhere 


(3) 


x( x ) 


0 

1 


In the domain f If the mask function vanishes identically and thus the original 
PDE in eq. (jT[) is satisfied. 

The difference between between the solution u of the original PDE and 
the solution of the penalized problem u^ (eq. ([ 2 J) ) is the penalization (or mod¬ 
eling) error, which depends on the size of the penalization parameter 77 , 

|| U-U V \\0QT] a (4) 

where a describes the order of the penalization method and || ■ || is a suitable 
norm, e.g., the energy norm || / || 2 = \f{x)\ 2 dx) l t' 2 or the L 1 or L°° norm. 

For convergence of the method the solution of the penalized problem u v should 
tend towards the solution of the original problem rt, i.e., lim^o || u — u v ||—» 
0. Thus a > 0 is required for convergence of the method. For the volume 
penalization we have a = 1/2 [lj[7]. 

Applying a numerical method to the penalized PDE eq. ([ 2 ]) we obtain the 
discretized penalized equation, 

L N u% = f N -^ X N (Bu-g) for x€D (5) 

with Aa; oc 1/N, where N denotes the number of grid points and L N is the 
discretized version of the operator L. In the case of an IBVP a suitable time 
discretization is also necessary. Explicit time discretization of the penalization 
term implies a time step restriction for stability. Typically the time step is 
limited by 77 and hence the problem becomes stiff in the limit 77 —>• 0. 

The solution of the discretized penalized equations , given that the dis¬ 
cretization is consistent and the numerical scheme is stable, converges towards 
the exact solution of the penalized equation (|2|. The discretization error, 

<8 

( 6 ) 

depends not only on the order of the numerical scheme, but is also limited by 
the regularity, i.e., the smoothness of the exact solution u v of the penalized 
problem. The order /3 is thus determined as the mimurn of the order of the 
underlying numerical scheme and the regularity of the exact penalized solution. 

Finally, the question in computations using the penalization method con¬ 
cerns the error between the exact solution u of the BVP (or IBVP) given by 
eq. 0 and the numerical solution of the penalized equation eq. ( |h| ) . This total 
error has two contributions: the modeling error which is due to the penaliza¬ 
tion; and the discretization error which is due to the numerical solution of the 
discretized problem. Applying the triangle inequality the following estimate 
holds: 

\\u-u% || < || u - uj| + || Ur, - u% || (7) 

and thus we obtain a bound for the total error. A straightforward argument 
would suggest choosing very small values for 77 to minimize the modeling error. 


u v - ||cx 
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However this would imply that the problem becomes stiff, one loses regularity in 
the exact solution of the penalized problem and the leading-order constant of the 
discretization error blows up and a very fine grid would become necessary. Thus 
to optimize the estimate both errors should be of the same order of magnitude. 
This shows that the penalization parameter r) and the numerical resolution N 
of the discretized problem are coupled and should be chosen accordingly. For 
further discussion we refer to [ 2 Uj . 



Figure 1: Sketch of the fluid domain f If immersed into the solid domain tt s . The 
boundary of the fluid domain is denoted by E = dflf. The total computational 
domain is S7 = U Q s . 


3 Volume penalization for Dirichlet and Neu¬ 
mann boundary conditions 

For illustration we consider first a simple toy problem, the Poisson equation in 
one space dimension, 

- = / for * e n f ( g ) 

completed with either homogeneous Dirichlet (u = 0 for x £ dVLf) or Neumann 
boundary conditions ( du/dx = 0 for x £ dflf). We choose Qf =]0, 7r[ and 
for the right-hand side / a trigonometric function, i.e., f(x) = m 2 sin mx and 
f(x ) = to 2 cosmx, respectively, with m £ N. The exact solution is u(x) = 
sin mx and u(x) = cos mx + C in the Dirichlet and Neumann case, respectively. 
The additive constant C £ R. shows that the solution is not unique for Neumann 
boundary conditions. To guarantee in this case the existence of a solution, the 
right-hand side / has to satisfy the compatibility condition, J Q f{x)dx = u'{x = 
7 r) — u' (x = 0) = 0. 

3.1 The penalized Poisson equation 

Now we replace the original problems by the penalized problems. The domain 
f =]0, 7r[ is embedded in the larger domain D =]0, 2-k[. Thus we have SI = 


5 



Clf U fl s , where f2 s is the penalization domain. For further simplification we 
impose periodic boundary conditions at the boundary of 0. 

In the Dirichlet case we obtain the penalized Poisson equation 

dj^u 1 

-7T+-^ = / fOT x e]0, 2 tt[ (9) 

ax A rj 

while in the Neumann case we have 

- ^ ((! - X') + VX) = / for x e ]°> 27T i ( 10 ) 

The penalization term in the Dirichlet case forces the solution to vanish inside 
the penalization domain and thus imposes homogeneous Dirichlet conditions at 
the interface. In the Neumann case the boundary condition corresponds to zero 
flux through the interface of fly and Q s . This can be achieved by imposing 
vanishing diffusivity (of order rj) inside the penalization domain. The function 
/ is extended in the larger domain by zero padding, i.e., f{x) = 0 for x £ fl s . 
Again both problems can be solved analytically in each sub-domain and we 
obtain in the Dirichlet case [291 


Ur,(x) 


sinmx + A\X + A 2 for x€]0,7t[ 

2 

x +^2 sin mx + Bi exp(-z/ jrj ) + B 2 exp(x/ y/rj) 


and in the Neumann case [T5] 


for 


x G]7r, 27 t[ 

( 11 ) 


J cos mx + A\X + A 2 for x e]0,7r[ 
\ BiX + B 2 for x g]7t, 27t[ 


( 12 ) 


Imposing continuity of the solution and of the derivative at x = 0 = 27t and 
x = 7T the coefficients (A 1; A 2 , Bi, B 2 ) can be determined. The corresponding 
values can be found in (25 1 , ITS] , Note that in the Neumann case only three out 
of the four coefficients can be determined as the solution is not unique. The 
penalization error || u — u v || can thus be explicitly computed and we find in 
the Dirichlet case the expected 0{y/rj) behavior, while in the Neumann case an 
0(rj) behavior is found which is better than O(yfrj) shown in ,T3| for the heat 
equation. In Fig. [ 2 ] (top) the exact solution of the penalized Dirichlet problem 
is plotted for two values of 77 (left) together with the first derivative (right). We 
find that for decreasing 77 the penalized solution tends towards the solution of the 
unpenalized problem. We also observe the existence of a boundary layer in the 
penalized domain (at the interface of flf and fl s , close to x = n) which becomes 
steeper when 77 becomes smaller. For 77 = 10~ 6 we find that the first derivative 
becomes almost discontinuous at x = tt. This shows that the regularity of u v is 
lost in the limit of small 77 and the problem becomes stiff. 


3.2 The discretized penalized Poisson equation 

To solve the penalized equations numerically we apply a second-order finite 
difference discretization to the penalized problems © and ( fl0| ). The compu¬ 
tational domain D = [0, 2 - 7 t ] is discretized with N grid points Xi = i/(2Tr),i = 
0,N — 1 applying periodic boundary conditions. We obtain the following 
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linear system in the Dirichlet case 


-D 2 + -xI)U = F 
V 


(13) 


with the vectors U = (u(x 0 ), ■ ■ ■, u(x jv-i)), X = (x( x o), ■ ■ •, xO*hv-i)), F = 
(f(x o)i • ■ •, f(x jv-i)) in R^, the identity matrix I and where 


D 2 


( ~ 2 1 

1 -2 

1 

h? 

V 1 


1 

1 


1 \ 

-2 1 
1 -2 j 


(14) 


is the second-order central finite difference operator with h = 2n/N. The re¬ 
sulting tridiagonal system can be solved using standard numerical linear algebra 
tools. 

In the Neumann case we get a second-order approximation with the following 
finite difference scheme, 


- * ( D f QD b + D b QD f ) U = F 


(15) 


where 0 = [9(x 0 ),9(x i),..., 0(ijv_i)] with 6{xi) = 1 — x( x i) + VX( x i))- The first 
order backward and forward finite difference operators Db and Dp are defined 
respectively, by 



( 1 



(- 1 1 

\ 

D b = \ 

-1 1 


, Dp = y 

-1 1 


h 


-1 1 ) 

h 

l 1 

-1 / 


(16) 

The linear system (15) is singular, the matrix has an eigenvalue 0 and a solution 
only exists if the right-hand side F is in its image. Special care has to be taken 
for solving the linear system using either the pseudo-inverse, or removing one 
equation. 

Performing numerical experiments |29l I18j showed that second-order con¬ 
vergence of the numerical solution of the penalized equations towards the ex¬ 
act solution can be obtained, under the condition that ry is sufficiently small. 
Moreover, for fixed ?y the total error (modeling plus discretization error) has 
a minimum at TV oc 1 /y/rj. Fig. [2] (bottom) illustrates this for the Dirichlet 
case. Varying ry with 1/TV like ry tx 1/TV 2 thus yields the best results in terms of 
minimum total error [29]. Nevertheless it must be mentioned that for spectral 
or fourth-order finite difference discretizations only first-order convergence can 
be observed. The explanation is that the linear system using second-order fi¬ 
nite differences is exactly the same for the penalized problem as for the original 
Dirichlet problem at all points inside the fluid domain, except at X\ and Xjv/ 2-1 
where values of O(ydy) are involved [29]. Hence for sufficiently small values of r] 
second-order convergence can be found. Note that for fourth order finite differ¬ 
ences the stencil is larger, thus more grid points are affected by the penalization 
term and the second-order convergence is lost. 
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Figure 2: Exact solution u v (top, left) and its first derivative u' v (top, right) of 
the penalized Poisson equation with Dirichlet boundary condtions for rj = ICE 2 
and 10 -6 . Convergence of the second-order finite difference scheme. Error with 
respect to the exact solution of the Dirichlet problem in f If (bottom). From [25] , 























3.3 The penalized Navier—Stokes and Maxwell equations 

Analogously to the penalized Poisson equation we can impose no slip and no 
penetration, i.e., Dirichlet boundary conditions in the incompressible Navier 
Stokes equations by adding a penalization term to the momentum equation, 

(9u 1 

^+uxu + vn - I/V 2 u = —x(u - u wa u) (17) 

at r\ 

V • u = 0 

where u = V x u is the vorticity and II = p + u 2 /2 is the modified pressure 
which satisfies — V 2 II = V • (uj x u+ -x( u ~ u wa u)). The kinematic viscosity is 
denoted by v. The velocity at the wall, u wa u does vanish identically in the case 
of fixed walls, or is prescribed in the case of moving walls. In the limit 77 —x 0 the 
solution of the penalized Navier-Stokes equations tends to the solution of the 
Navier-Stokes equations with no-slip boundary conditions and the penalization 
error is of order 0(y/rj), as shown in [I, 7j. 

The induction equation which describes the evolution of the magnetic field 
B can be penalized in a similar way and we obtain 

1 

— - V x (u x B) - AV 2 B = —x(B - B wall ) (18) 

at ij 

V • B = 0 

where A is the magnetic diffusivity. The magnetic field at the wall, B wa u , can be 
freely chosen and hence not all components have to be penalized. For example 
choosing B wa u = B\\ where B\\ is the component of B parallel to the wall, only 
penalizes the normal component and leaves the parallel component free. This 
allows modeling perfectly conducting boundary conditions. Note that in the 
MHD case no mathematically rigorous convergence theorem has been proven 
yet, but in | 2 fij asymptotic arguments for estimating the penalization error have 
been given. 

For the transport of a passive scalar £ we consider the advection-diffusion 
equation where no-flux conditions V£ • n = 0, i.e., homogeneous Neumann 
boundary conditions are imposed [13], 

d£ 

+ [(1 - X)u] • V£ = V • [«(1 - x) + (19) 

where 77 ^ is the penalization parameter for the scalar £ and k its diffusivity. In 
the limit 77 c —x 0 the solution of the penalized equation tends to the solution 
of the advection-diffusion equation with no-flux boundary conditions and the 
penalization error is of order 0(y/rj), as shown in [13| . 

4 Application to fluid turbulence 

In the following we consider different applications to two- and three dimensional 
hydrodynamic incompressible flows with no-slip walls and also the transport of 
a passive scalar where no-flux boundary conditions are imposed. The numerical 
method is based on a pseudo-spectral discretization of the penalized Navier 
Stokes and advection-diffusion equations and for details we refer the reader 

to [Ml nuns EH- 
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4.1 Two-dimensional turbulence in a circular container 

In [39j we presented numerical simulations of two-dimensional decaying tur¬ 
bulence in a circular container with no-slip boundary conditions computed at 
resolution 1024 2 using the volume penalization technique. Starting with random 
initial conditions with Reynolds number Re = 5 x 10 4 where Re is based on 
the domain size, the turbulent kinetic energy and the kinematic viscosity v, the 
flow rapidly exhibits self-organization into coherent vortices. Two snapshots of 
the vorticity field oj at later times are shown in Fig. [3] One-dimensional cuts 
in Fig. [4] (top) illustrate the intermittent character of the vorticity field and 
the spikes at the domain boundary show the strong production of vorticity in 
a thin boundary layer due to the no-slip boundary conditions. The cut of the 
mask function together with the cuts of the velocity components confirm that 
in the penalization domain the velocity does indeed vanish. The formation of 
coherent vortices and the viscous boundary layer have significant impact on the 
production and decay of integral quantities. The evolution of kinetic energy 
E{t), enstrophy Z(t) and palinstrophy P(t) are shown in Fig. [4] (bottom). The 
corresponding balance equations read, 

d t E = — 2uZ 1 d t Z = —2vP + ® u> (n • Vw) ds , (20) 

Jan 

where n denotes the outer normal with respect to dfl. The source term in 
the enstrophy dissipation equation involves the vorticity and its gradient at 
the boundary and yields a significant contribution in the small viscosity limit. 
The no-slip wall produces vortices which are injected into the bulk flow and 
tend to compensate the enstrophy dissipation as observed in Fig. [4] The self- 



Figure 3: Decaying two-dimensional flow in a container with no-slip boundary 
conditions at initial Reynolds number 5 x 10 4 . Snapshots of vorticity w. From 

m- 


organization of the flow is also reflected by the transition of the initially Gaussian 
vorticity probability density function (PDF) towards a distribution with expo¬ 
nential tails. Because of the presence of coherent vortices the pressure PDF 
becomes strongly skewed with exponential tails for negative values. Details can 
be found in [55] . 
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Figure 4: Decaying two-dimensional container with no-slip boundary conditions 
at initial Reynolds number 5 x 10 4 . Id vertical cuts of vorticity, the two velocity 
components and the mask function (top). Evolution of the kinetic energy E{t) = 
\ J \u\ 2 dx : enstrophy Z(t) = | f \u\ 2 dx and palinstrophy P(t ) = \ f \S7u\ 2 dx 
in double logarithmic representation (bottom), where r = t/t e is based on 
the initial eddy turn over time t e = yj2Z(t + 0) = 0.061. The flow has been 
integrated for 650t e , corresponding to more than 10 5 time steps. The inset 
shows the corresponding log-lin representation. From }39| . 
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4.2 Passive scalar transport in two-dimensional confined 
turbulence 


To illustrate the volume penalization for imposing no-flux boundary conditions 
we show a numerical simulation of a flow with a passive scalar in a simplified 
mixing device [13]. The penalized Navier-Stokes equations (18) are solved to¬ 
gether with the penalized advection-diffusion equation (19) in a square domain 
with periodic boundary conditions using a Fourier pseudo-spectral method. The 
fluid domain corresponds to a circular vessel in which a cross-shaped rotor in the 
center of the domain rotates in the clockwise direction. The boundary condi¬ 
tions are no slip for the velocity and no flux for the passive scalar. The vorticity 



Figure 5: Vorticity field in a circular vessel with a cross-shaped clockwise 
rotating rotor (left). Corresponding passive scalar field (right). From [15] . 

field in Fig. [5] illustrates the formation of boundary layers at the wall and at 
the rotor. The formed vortex sheets destabilize and detach forming coherent 
vortices which are ejected into the bulk flow. The corresponding passive scalar 
field, the initial condition corresponds to a Gaussian blob, is advected by the 
mean rotation induced by the rotor. The mixing process is further enhanced by 
the coherent vortices generated by the rotor blades and the boundary layer of 
the vessel. For further details we refer the reader to [T5] , 

4.3 Flow past flexible flapping plates in three dimensions 

Motivated by simplified models for swimming organisms or robots, which rely 
on chord-wise flexible elastic plates, we present numerical simulations of fluid- 
structure interaction using the volume penalization. We consider a plate made 
out of linearly elastic inextensible material, which is perfectly rigid in the span- 
wise direction but flexible in the chordwise direction. The plate can be thus 
modeled by the nonlinear beam equation with clamped free boundary condi¬ 
tions and it is solved with classical finite differences. The fluid part is solved 
with a pseudo-spectral method and volume penalization m • Depending on 
the fluid/plate density ratio up to 25 iterations of the fluid-solid coupling are 
necessary within each time step. Details on the numerical method are described 
in EH- ' 
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In [101 we considered a configuration with imposed mean flow, and imposed 
at the leading edge of the plate a sinusoidal pitching motion. The Reynolds num¬ 
ber is about Re « 1000. We first simulated a swimmer with a rectangular plate 
and compared the results with a recent experimental study, before considering 
also an expanding and a contracting shape of the plate. Flow visualizations for 
the three geometries showing vorticity isosurfaces are given in Fig. [6] The tip 
vortices observed in Fig. [6] originate from three-dimensional effects due to the 
finite span. These have important effects for predicting the swimmers cruising 
velocity, since they contribute significantly to the drag force. We found that the 
cruising velocity of the contracting swimmer is larger than of the rectangular 
one, which in turn is larger than the expanding one. This observation can be ex¬ 
plained by the relative importance of the tip vortices which interact differently 
with the flexible plates for the three considered geometries of the swimmer. For 
the contracting case the tip vortices rapidly detach and thus reduce drag, while 
in the expanding case they are attached down to the trailing edge. 



Figure 6: Flow generated by flapping cord-wise flexible plates of contracting 
(left), rectangular (center) and expanding shape (right) with imposed pitching 
motion at their leading edge and imposed mean flow (from left to right). Shown 
are isosurfaces of vorticity || w ||= 17.5. From m- 


5 Applications to plasma turbulence 


For the numerical simulation of plasma turbulence we consider the non-ideal 


MHD equations (eqs. 18 and 19) in which both viscous and resistive effects are 
taken into account. The magnetic Prandtl number, defined as the ratio between 
kinematic viscosity v and magnetic diffusivity A, is equal to one. An isothermal, 
incompressible plasma is considered with uniform and constant transport coef¬ 
ficients. This approximation simplfies the problem as much as possible, while 
retaining the required level of complexity to study the nonlinear dynamics. The 
boundary conditions corresponding to solid domains which are perfect conduc¬ 
tors are imposed with the volume penalization method. The numerical code is 
based on a Fourier pseudo-spectral discretization using FFTs to compute the 
derivatives and to solve the Poisson equations. It is described in detail in [25] . 
including benchmarking and detailed validation studies. 
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5.1 Spin-up in two-dimensional confined MHD 

We consider first two-dimensional decaying MHD turbulence in bounded do¬ 
mains and investigate the spontaneous self-organization with a particular em¬ 
phasis on the symmetry-breaking induced by the shape of the confining bound¬ 
aries; for details we refer to [3] and m- This symmetry-breaking is quantified 
by the angular momentum, which is shown to be generated rapidly and sponta¬ 
neously from initial conditions free from angular momentum when the geome¬ 
try lacks axisymmetry. In [5] this effect was illustrated by considering circular, 
square, and elliptical boundaries. It was shown that the generation of angular 
momentum in non-axisymmetric geometries can be enhanced by increasing the 
magnetic pressure. Moreover, the effect becomes stronger at higher Reynolds 
numbers, which are based on the turbulent kinetic energy, the length scale of 
the domain and the kinematic viscosity. The generation of magnetic angular 
momentum or angular field, previously observed in j4] at low Reynolds numbers, 
becomes weaker at larger Reynolds numbers. For axisymmetric geometries, the 
generation of angular momentum is absent; nevertheless, a weak magnetic field 
can be observed. The derived evolution equations for both the angular momen¬ 
tum and angular field yield possible explanations for the observed behavior. 



Figure 7: Two-dimensional magnetized plasma: stream function (left) and 
vorticity (right) for the ovoid (top) and the D-sliaped geometry (bottom). The 
visualizations correspond to the time-instants at which the absolute value of the 
angular momentum reaches its maximum. From Ell¬ 
in Fig. [7] we show simulations of two-dimensional decaying MHD turbu¬ 
lence inside an ovoid and in a D-shaped geometry starting with random initial 
conditions ED- For both geometries we observe that the two-dimensional mag¬ 
netized plasma self-organizes into a state containing large-scale flow structures, 
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time 


Figure 8: Two-dimensional magnetized plasma: Time evolution of the angular 
momentum L u (t) for the ovoid and the D-shaped geometry. From 141] . 


illustrated by the stream function in Fig. [7] (left) and vorticity, Fig. 0 (right). 
To quantify the spin-up we consider the angular momentum defined as 

L u (t) = / e z ■ r x uds = 2 / i/j ds (21) 

L u quantifies the fluid rotation. The maximum angular momentum for a given 
kinetic energy is obtained for a fluid in solid body rotation. In this particular 
realization we find that the generation of angular momentum is stronger in the 
ovoid than in the D-shaped geometry, as shown in Fig. [8] 


5.2 Self-organization of confined MHD flows in toroidal 
domains 

The spatio-temporal self-organization of visco-resistive magnetohydrodynamics 
in a toroidal geometry (see Fig. 10 left) was studied in [24] using fully three- 
dimensional simulations considering two geometries: a torus with a symmetric 
poloidal cross-section and one with an asymmetric poloidal cross-section. The 
magnetized plasma is initially in a quiescent state and curl-free toroidal magnetic 
and electric fields are imposed. The simulations show that spontaneously a 
flow field is generated in both geometries and the magnetized plasma starts to 
move. Moreover, the flow evolves from dominantly poloidal to toroidal when the 
Lundquist (or Reynolds) numbers M are increased. Here the viscous Lundquist 
number M is defined as M = Ca^/v where Ca is the toroidal Alfven velocity, 
L the diameter of the cross section of the torus and v the kinematic viscosity. In 
[24] we have shown that this toroidal organization of the flow is consistent with 
the tendency of the velocity field to align with the magnetic field. Furthermore, 
we found that the up-down asymmetry of the geometry causes the generation 
of a non-zero toroidal angular momentum. 

Figure [9] illustrates the streamlines, colored with the toroidal velocity value 
for two Lundquist numbers for the asymmetric poloidal geometry. For the low 
Lundquist number case (Fig. [9j left) we do indeed observe a pair of counter¬ 
rotating vortices in the poloidal plane, while for the larger Lundquist number 
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Figure 9: Confined MHD flows in toroidal domains. Asymmetric poloidal 
geometry. Streamlines colored with toroidal velocity (uq) for M = 7.5 (left) and 
M = 75.2 (right). From [24]. 




(Fig.|9j right) the flow starts moving in the toroidal direction. A similar behavior 
is observed for the symmetric poloidal cross-section. 



Figure 10: (a) Sketch of the toroidal geometry with asymmetric poloidal cross- 
section, (b) toroidal velocity field component isocontours (blue +9 • 10 4 , orange 
9 • 10 4 ) and (c) perturbed toroidal magnetic field component isocontours (red 
+0.025, orange +0.04, yellow +0.05). From [25]. 


The velocity and the perturbed magnetic toroidal component isocontours at 
the steady state for a toroidal geometry with asymmetric poloidal cross section 
are shown in Fig. 10 (b and c). The perturbed toroidal magnetic field is created 
by the velocities in the poloidal plane. This component of the magnetic field is 
important because it generates a toroidal Lorentz force that induces the toroidal 
velocities. 

Figure 11 (left) shows that the toroidal velocity increases with the viscous 
Lundquist number in both geometries and saturates at about 86 % of the total 
squared speed. The inset shows the modulus of the cosine of the angle be¬ 
tween the velocity and the magnetic field and thus quantifies that the velocity 
fields tend to align with the magnetic field for increasing M. However, for the 
volume-averaged toroidal angular momentum defined as ( Lg) = y f v Rug dV, 
we observe fundamental differences for the two geometries as shown in Fig. ITT] 
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(right). For the torus with the symmetric poloidal cross-section (Lg) identi¬ 
cally vanishes for all considered Lundquist numbers. In contrast we observe for 
the asymmetric case that the toiroidal angular momentum increases with the 
Lundquist number. 



M 
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Figure 11: Confined MHD flows in toroidal domains. Left: normalized toroidal 
velocity (ug) / ( u 2 ). The inset shows the modulus of the cosine of the angle 
between the velocity and the magnetic field. Right: Volume averaged toroidal 
angular momentum (Lg) . All quantities are plotted as a function of the viscous 
Lundquist number M for the symmetric and asymmetric poloidal cross section 
geometries. From [21j . 


5.3 Effect of toroidicity in RFP dynamics 

Finally, we consider the reversed field pinch (RFP) dynamics and study the in¬ 
fluence of the curvature of the imposed magnetic field. In RFP experiments the 
plasma evolved to quasi-stationary equilibria characterized by Beltrami mini¬ 
mum energy states for which the magnetic field corresponds to eigenfunctions 
of the curl operator [30]. In SB] we compared the RFP flow of a magnetofluid 
in a torus with aspect ratio 1.83 with the flow in a periodic cylinder and studied 
the persistence of these dominant helical modes for varying pinch ratios. The 
flow configuration of both geometries is illustrated in Fig. [T2] The pinch ratio 
is defined as the wall-averaged poloidal magnetic field divided by the volume- 
averaged toroidal magnetic field, 0 = Bp/(Bt) ■ The ratios of kinetic energy 
of the dominant axial and toroidal mode of the total kinetic energy versus the 


pinch ratio are shown in Fig. 13 (top) and (bottom) for the cylindrical and the 
toroidal geometry, respectively. We find that an axisymmetric toroidal mode 
is always present in the toroidal, but absent in the cylindrical configuration. 
In particular, in contrast to the cylinder, the toroidal case presents a double 
poloidal recirculation cell with a shear localized at the plasma edge. Quasi- 
single-helicity states are found to be more persistent in toroidal than in periodic 
cylinder geometry. The dominant helical modes at pinch ratios 0 > 2 are illus¬ 
trated in Fig. [13] by showing axial and toroidal velocity isosurfaces to get further 
insight into the flow topology. For further details we refer to [2BI . 
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Figure 12: Flow configuration of the RFP dynamics. Left: toroidal geometry. 
Right: periodic cylinder. From |2(jj . 

6 Conclusion and perspectives 

We reviewed immersed boundary methods with a special focus on the volume pe¬ 
nalization method for imposing Dirichlet (corresponding to no-penetration and 
no-slip conditions) or Neumann boundary conditions (corresponding to no flux 
conditions) in complex geometries. The mathematical concepts for choosing the 
parameters involved, i.e., the penalization parameter and the grid size have been 
illustrated considering simple one dimensional problems. Applications to hydro- 
dynamic and magnetohydrodynamic turbulence in complex geometries using a 
simple Fourier pseudo-spectral method, which can be parallelized on massively 
parallel machines using standard libraries like P3DFFT, illustrated the versatile 
use of this technique for various problems encountered in computational physics. 
Toroidal geometries can thus be efficiently handled even including asymmetric 
poloidal cross-sections and simulations for higher Lundquist numbers become 
feasible. An essential feature of the volume penalization method is that it be¬ 
comes more attractive for computing fluid flows for small viscosity values, i.e., 
for high Reynolds/Lundquist number flows. The reason is that the effective 
penalization boundary layer size depends on the product of viscosity and the 
penalization parameter and thus the method becomes more precize without 
using prohibitively small penalization parameters, cf. |29| . 

One perspective of current research is the development of higher-order penal¬ 
ization methods which allow faster convergence to be obtained. The Cartesian 
grid introduces in dimension larger equal to two a staircase effect for complex 
(non grid aligned) geometries and the approximation of the mask function thus 
reduces to first order. Techniques based on interpolation to obtain higher order 
for complex geometries have been proposed, e.g., in [5B] in the context of finite 
volume formulations. Another challenging topic are more sophisticated bound¬ 
ary conditions for MHD flows overcoming the limitation of perfect conductors 
in the surrounding solid domain. 


18 












Figure 13: RFP dynamics. Ratio of the kinetic energy of the dominant axial 
(top)/toroidal (bottom) modes over the total kinetic energy for the cylindrical 
(top) and the torus geometry (bottom) for M = 888 as a function of the pinch 
ratio 0. Visualization of the modes: Top: axial velocity isosurfaces +0.008 
(blue) and 0.008 (orange). Bottom: toroidal velocity isosurfaces +0.007 (blue) 
and 0.007 (orange). From [26] . 
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